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Abstract.  The  problem  of  the  consolidation  of  an  aerated  fine  powder  under  gravity  is  consid¬ 
ered.  The  industrial  relevance  of  the  problem  is  discussed  and  a  mathematical  model  is  introduced. 
The  mathematical  structure  is  that  of  a  coupled  system  for  three  unknowns,  pressure,  stress  and 
height  of  the  powder  in  the  (axisymmetric)  bunker  containing  it.  The  system  itself  consists  of  a 
parabolic  PDE,  an  ODE  and  an  integral  equation  determining  a  free  boundary  corresponding  to  the 
height  of  the  powder.  Existence  and  uniqueness  of  a  solution  is  established.  A  numerical  method 
based  on  a  formulation  of  the  semidiscretized  problem  as  an  index  1  DAE  is  proposed  and  imple¬ 
mented.  The  feasilibility  of  the  approach  is  illustrated  by  computational  results. 
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1.  Introduction.  One  important  factor  determining  the  mechanical  properties 
of  fine  powders  is  the  possible  presence  of  an  interstitial  fluid,  say  air.  Indeed,  any 
pressure  gradient  clearly  translates  into  an  additional  body  force.  In  gravity  flow, 
this  force  may  be  of  the  order  of  the  weight  density.  One  then  speaks  of  fluidization, 
and  the  air-powder  system  essentially  behaves  as  a  liquid.  This  paper  is  concerned 
with  the  “opposite”  phenomenon:  consolidation.  More  precisely,  consider  storing 
some  fine  powder  in  a  bunker  or  silo,  see  Figure  1.1.  Inevitably,  some  air  will  get 
trapped  during  filling.  The  corresponding  partial  fluidization  can  have  very  serious 
and  unwanted  consequences  in  practice  and  may  result,  upon  retrieval  of  the  material, 
in  uncontrollable  flows  and  flooding.  The  excess  air  does  diffuse  through  the  powder 
and  eventually  escapes  through  the  top  surface.  A  natural  question  is  then:  how  long 
does  one  have  to  wait  for  the  air-powder  system  to  settle  and  allow  for  safe  handling? 
This  is  the  motivation  of  this  paper.  We  note  that  similar  questions  may  apply  in 
Soil  Mechanics  in  general  and  the  study  of  landfills  in  particular. 

By  a  powder,  we  mean  a  material  consisting  of  many  individual  solid  particles  of 
sizes  roughly  between  10~7m  and  10~5m  [10].  If  the  particles  are  much  larger,  then 
the  gas  can  circulate  nearly  freely  between  them  and  thus  is  not  an  important  factor. 
Although  the  present  study  could  be  easily  generalized  to  any  type  of  axisymmetric 
container,  we  consider  for  the  sake  of  simplicity  a  vertical  cylindrical  bunker  containing 
a  granular  material  subject  to  gravity,  see  again  Figure  1.1.  We  work  under  several, 
more  restrictive,  simplifying  assumptions.  First,  the  problem  is  made  essentially  one¬ 
dimensional  by  assuming  all  the  physical  variables  to  be  uniform  across  horizontal 
sections.  Such  a  simplification  is  clearly  not  fully  justifiable  in  general.  However,  it 
has  been  shown  to  lead  to  meaningful  asymptotic  results  in  the  context  of  a  so-called 


‘Department  of  Mathematics  and  Center  for  Research  in  Scientific  Computation,  North  Car¬ 
olina  State  University,  Raleigh,  NC  27695-8205,  USA  (gremaud@math.ncsu.edu).  Partially  supported 
by  the  Army  Research  Office  (ARO)  through  grants  DAAH04-95-1-0419,  DAAH04-96-1-0097  and 
DAAD19-99-1-0188,  by  the  National  Science  Foundation  (NSF)  through  grant  DMS-9818900,  and 
by  a  grant  from  the  North  Carolina  Supercomputing  Center. 

^Department  of  Mathematics  and  Center  for  Research  in  Scientific  Computation,  North  Car¬ 
olina  State  University,  Raleigh,  NC  27695-8205,  USA  (ctk@math.ncsu.edu).  Partially  supported  by 
National  Science  Foundation  grant  DMS-9700569. 

*Jenike  &  Johanson,  inc,  Westford,  MA  01886,  USA  (royalta@jenike.com). 

^Department  of  Mathematics  and  Center  for  Research  in  Scientific  Computation,  North  Carolina 
State  University,  Raleigh,  NC  27695-8205,  USA  (kahyman@eos.ncsu.edu).  Partially  supported  by  a 
Department  of  Education  GAANN  fellowship. 


1 


Report  Documentation  Page 

Form  Approved 

OMB  No.  0704-0188 

Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 

VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 

1.  REPORT  DATE 

2000 

2.  REPORT  TYPE 

3.  DATES  COVERED 

00-00-2000  to  00-00-2000 

4.  TITLE  AND  SUBTITLE 

On  a  Powder  Consolidation  Problem 

5a.  CONTRACT  NUMBER 

5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

6.  AUTHOR(S) 

5d.  PROJECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

North  Carolina  State  University, Center  for  Research  in  Scientific 
Computation, Raleigh, NC, 27695-8205 

8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

10.  SPONSOR/MONITOR'S  ACRONYM(S) 

11.  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 

12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited 

13.  SUPPLEMENTARY  NOTES 

The  original  document  contains  color  images. 

14.  ABSTRACT 

see  report 

15.  SUBJECT  TERMS 

16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION  OF 

ABSTRACT 

18.  NUMBER 

OF  PAGES 

19 

19a.  NAME  OF 

RESPONSIBLE  PERSON 

a.  REPORT 

unclassified 

b.  ABSTRACT 

unclassified 

c.  THIS  PAGE 

unclassified 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std  Z39-18 


2 


P.A.  GREMAUD,  C.T.  KELLEY,  T.A.  ROYAL  AND  K.A.  COFFEY 


Fig.  1.1.  Geometry  and  coordinate  systems  for  the  vertical  cylindrical  bunker.  The  height  of 
the  column  of  powder  of  time  t  is  denoted  H(t). 


Janssen  analysis  [6]  of  the  behavior  of  columns  of  granular  materials,  [10],  p.84-90. 
The  present  study  can  in  fact  be  viewed  as  a  generalization  of  Janssen’s  original 
approach  to  fine  powders,  where  the  presence  of  air  cannot  be  neglected. 

The  model  that  is  derived  in  Section  2  directly  results  from  classical  conservation 
principles  together  with  Darcy’s  law.  The  unknowns  are  the  pressure  of  the  gas  p, 
the  vertical  stress  a  and  the  height  of  the  column  of  powder  H.  The  mathematical 
structure  of  the  problem  is  nonstandard  as  it  consists  of  a  system  of  three  equations: 
a  parabolic  PDE,  an  ODE  and  an  integral  equation,  which  corresponds  to  a  nonlocal 
equation  for  the  top  free  boundary.  An  earlier  and  slightly  different  model  was  derived 
in  [1],  In  Section  3,  a  much  simpler  auxiliary  “toy  problem” ,  that  roughly  corresponds 
to  heat  conduction  in  an  expanding  rod  is  considered,  see  also  [9]  for  a  numerical  study 
of  a  closely  related  case.  That  problem  shares  some,  but  not  all,  of  the  difficulties  of 
the  full  problem,  and  its  analysis  is  meant  to  provide  a  road  map  for  what  is  covered 
in  the  rest  of  the  paper.  Section  4  is  devoted  to  the  analysis  and  numerical  analysis 
of  the  full  problem.  Existence  and  uniqueness  of  a  solution  is  established.  A  simple 
numerical  method  is  proposed.  Essentially,  the  height  H  is  expressed  as  a  function 
of  the  other  variables.  This  expression  is  used  to  transformed  the  problem  into  one 
defined  in  a  fixed  spatio-temporal  domain.  Then,  the  remaining  system,  involving 
now  two  unknowns,  is  semidiscretized  in  space.  This  results  in  a  DAE  which  is  solved 
through  the  use  of  a  linearized  Implicit  Euler  method.  The  feasibility  and  efficiency  of 
the  method  is  illustrated  in  Section  5  where  computational  experiments  are  presented. 
Finally,  some  closing  remarks  are  offered  in  Section  6. 


2.  The  model.  We  denote  by  T  the  solid  density,  i.e.,  the  actual  density  of 
the  particle;  T  is  assumed  to  be  constant.  The  density  p  of  the  gas  is  an  unknown 
function  of  the  time  t  and  the  position  z,  i.e.  the  height,  see  Figure  1.1),  under  the 
above  simplifying  assumption.  A  most  important  quantity  is  the  bulk  density  denoted 
by  7.  The  bulk  density  can  be  defined  as  the  effective  density  of  the  granular  material. 
More  precisely,  if  fs  stands  for  the  volume  fraction  occupied  by  the  solid,  we  have 


(2.1) 


7  =  fBT  +  (1  -  f„)p. 
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Fig.  2.1.  Graph  of  the  bulk  density  7  as  a  function  of  the  stress  a . 


Typically,  the  gas  density  p  is  at  least  three  orders  of  magnitude  less  than  T.  Conse¬ 
quently,  the  solid  fraction  fs  can  in  fact  by  approximated  by 


7 

r' 


The  bulk  density  is  found  to  increase  with  the  application  of  increasing  stresses.  The 
following  relation  is  well  supported  by  experimental  evidence  [7] 


(2-2)  7  =  7™(1  +  — 

where  0  <  f3  <  1  is  the  coefficient  of  compressibility,  7m  >  0  and  am  >  0  being  two 
material  constants  and  where  a  stands  for  the  vertical  stress.  Note  that  above  formula 
makes  sense  only  for  “reasonable”  values  a.  Indeed,  one  does  not  expect  7  to  increase 
without  bounds  for  increasingly  large  stresses.  For  instance  7  should  certainly  not 
exceed  T.  Accordingly,  the  mathematical  analysis  of  the  problem  was  be  done  under 
the  following  assumptions  on  7,  see  Figure  2.1,  which  are  consistent  with  (2.2)  for  a 
not  too  large 


(2.3)7  e  C°°{R),  7(0)  =  7m  >  0,  0  <  7(c)  <  7  <  r,  0<7,(a)<7(x),  Va. 


In  what  follows,  the  cylindrical  container  is  assumed  to  be  of  constant  circular  cross- 
section  A,  but  this  is  not  essential,  see  §6.  The  spatial  coordinates  are  taken  as  the 
cylindrical  coordinates  r  and  z,  see  Figure  1.1.  Under  the  assumption  of  horizontal 
uniformity,  no  angular  variable  is  needed.  We  denote  by  H[t)  the  height  of  the  column 
of  powder  at  time  t.  Since  no  powder  escapes  during  the  consolidation  process,  the 
total  mass  M  of  solid  in  the  bunker  is  conserved.  Under  the  above  approximation  of 
fs,  the  global  conservation  of  solid  reads 


(2-4) 


t  >  0, 


where  R  denotes  the  radius  of  the  container.  The  above  equation  can  be  considered 
as  “the”  equation  for  H,  i.e.,  the  relation  defining  the  top  free-boundary.  Its  non-local 
character  should  be  noted.  Local  conservation  of  solid  and  gas  yield 


dtl  +  dz(jus)  =  0, 
9t{{  1  -  ^)p)+dz(p(  1  -  ^)ug)  =  0, 


(2.5) 

(2.6) 
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where  us  and  ug  stand  for  the  velocity  of  solid  and  gas  respectively. 

The  gas  is  assumed  to  be  ideal  and  isothermal.  In  other  words,  if  p  denotes  the 
pressure  of  the  gas,  we  have 


(2.7) 


P  =  Po 
P  Po’ 


where  po  and  po  are  two  constant  reference  values.  Velocities  and  pressure  are  related 
to  each  other  through  Darcy’s  law,  i.e., 


(2.8)  ug-us  =  -K(^)dzp, 

where  K  =  K( 7)  is  the  permeability.  The  gas  flow  resulting  from  a  pressure  gradient 
is  obviously  also  dependent  on  the  bulk  density.  Here,  we  take 


(2.9)  K( 7)  =  K0  (-l)  , 

where  Kq  and  70  are  reference  values  and  a  is  a  positive  constant.  The  parameters  /3, 
am,  7 m,  a,  Ko  and  70  appearing  in  (2.2)  and  (2.9)  can  be  determined  experimentally. 

Given  the  symmetry  of  the  problem  as  considered,  the  stress  tensor  T  has  the 
form 


T  = 


(T'p'p  ^  v  z 
®rz  O zz 


In  first  approximation,  we  will  assume  that  the  vertical  and  horizontal  stresses  are 
principal  stresses.  This  amounts  to  neglecting  the  contribution  from  arz  in  T,  again 
a  questionable  assumption  in  some  cases  [10],  p.84-90. 

Next,  consider  a  balance  of  forces  acting  on  an  infinitesimal  slice  of  radius  R  and 
height  5z  of  the  material,  see  again  Figure  1.1.  The  various  forces  are 

•  weight  of  solid:  —jn  R2  5z\ 

•  if  tw  is  the  wall  shear  stress,  there  is  an  upward  force  of  2  7r  R  5z  tw ; 

•  pressure  at  bottom:  p(z),  and  top:  —  (p(z)+6p);  this  creates  a  force  —ttR2  Sp; 

•  stress  at  bottom:  azz(z),  and  top:  —(azz(z)  +  8azz);  (compressive  stresses 
are  taken  as  positive  for  granular  material);  this  creates  a  force  —ttR2  Sazz. 

The  resulting  force  balance  gives 


dzazz  +  dzp  -  —tw  +7  =  0. 

Further,  by  applying  the  law  of  sliding  friction  on  the  wall,  one  finds 

T~w  —  PwCrr, 

where  pw  is  the  coefficient  of  wall  friction.  Finally,  the  two  remaining  components  of 
the  stress  tensor  T  are  related  through  a  plasticity  model.  The  powder  is  taken  to  be 
an  ideal,  cohesionless,  Coulomb  material  [10].  Since  we  assume  vertical  and  horizontal 
stresses  to  be  principal  stresses,  this  implies  that  their  ratio  has  to  be  constant.  More 
precisely,  we  have 


arr  1  sin  S 

azz  1  ±  sin  5 


(2.10) 
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where  6  is  the  angle  of  internal  friction,  i.e.,  the  coefficient  of  internal  friction  p  is 
given  by  p  =  tan  5.  The  angle  5  is  also  equal  to  the  angle  of  repose  of  the  material. 
The  plus  or  minus  signs  in  (2.10)  correspond  to  either  the  realization  of  the  active 
state  if  J  =  ,  i.e.,  azz  >  arT ,  or  the  passive  state  if  J  =  ,  i.e.,  <jzz  <  <jrr. 

We  will  work  here  under  the  assumption  that  the  material  is  in  the  active  state,  since 
this  state  is  the  one  observed  upon  filling.  This  however  has  very  little  consequence 
on  the  rest  of  the  study  (change  of  the  value  of  the  constant  J) . 

Using  both  (2.10)  and  the  form  of  tw  given  above,  we  can  eliminate  o>r  from  the 
equations.  Denoting  by  a  the  remaining  stress  component  crZ2,  we  get 

(2.11)  dza  +  dzp-Kcr+  7  =  0, 

where  the  constant  k  is  defined  by  k  =  2pw  J/R.  We  now  eliminate  the  velocities 
from  the  system.  Relations  (2.5)  and  (2.8)  combine  into 

da  +  dz  (7  ug  +  7  Kdzp)  =  0. 

After  integration  in  space,  and  taking  into  account  the  boundary  conditions  ug( 0,  t)  = 
0  and  dzp(0,  t)  =  0,  we  obtain 


1  fz 

ug  =  -~  da  ~  K dzp. 

7  Jo 

Next,  using  (2.7),  one  can  rewrite  (2.6)  in  terms  of  the  pressure  p,  eliminating  p  from 
the  problem.  Plugging  the  expression  for  ug  in  (2.6)  then  yields  our  last  missing 
equation 

(i-l)dtP-^da-dz  (p(^-f))  JqZ  da-(i-^)dz(pKdzp)  +  ^dadzp  =  o. 

One  last  simplification  is  considered.  The  last  term  in  the  previous  relation  has  very 
little  influence  on  the  problem.  Although  not  obvious  analytically,  this  was  carefully 
verified  numerically:  for  realistic  values  of  the  various  parameters,  the  solution  changes 
by  less  than  .1%  when  switching  this  term  on  and  off.  Omitting  this  term  simplifies 
somewhat  the  analysis  on  the  one  hand,  and  does  not  appear  to  influence  the  solution 
in  any  observable  way  on  the  other  hand.  In  conclusion,  and  under  the  previous 
remarks  and  assumptions,  the  unknowns  a,  p  and  H  are  to  be  determined  by  the 
following  system  of  three  integrodifferential  equations  in  {(z,t);  0  <  t,  0  <  z  <  H(t)} 


(2.12) 

(2.13) 

(2.14) 


^)dz(pKdzp)  =  0, 


Those  equations  are  combined  with  the  initial  and  boundary  conditions 


(2.15)  pM)=po, 

(2.16)  a(H(t),t)  =  0,  <9zp(0,f)=0,  p(H(t),t)  =  patm, 


where  patm  stands  for  the  value  of  the  atmospheric  pressure. 
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3.  An  auxiliary  problem.  To  guide  us,  as  well  as  the  reader,  through  the  anal¬ 
ysis  of  the  above  problem,  an  auxiliary  problem  that  shares  some  of  the  nonstandard 
features  encountered  in  the  previous  section  is  first  analyzed.  The  problem  below  can 
be  regarded  as  the  description  of  the  behavior  of  a  one-dimensional  rod  subject  to 
thermal  expansion.  Although  such  an  interpretation  is  not  strictly  needed,  one  can 
think  of  the  unknowns  appearing  below  as  the  temperature  9,  the  linear  density  of 
the  rod  p  and  its  length  at  time  t,  s(t).  The  problem  is  then 


rs(t) 

Jo 


dt0  -  dxx9  =  0 
9(0,  t)  =  9(s(t),  t)  =0 

9(x,  0)  =  90(x) 

p(9(x,  t))  dx  =  M  =  1 


0  <  x  <  s(t),  t  >  0 
t  >  0, 

0  <  x  <  s(0), 
t  >  0, 


where  9o  is  a  given  initial  condition.  We  start  by  mapping  the  problem  into  a  fixed 
spatial  domain  by  introducing  the  variables 


y  =  x/s(t)  u(y,  t)  =  9(x,  t). 
The  heat  equation  then  turns  into 

ys'(t)  1 

9tU  S(t)  dyU  s(tydvvU~Q- 
Now,  f0s(t)  p(9(x,  t))  dx  =  fg  p(u(y,  t ))  s(t)dy  =  1  and  thus 


* (*)  =  7i — - — —  311(1  «'W 

Jo  P(u(yi  t))  dy 


Finally,  we  obtain 


~  fo  P' (u)dtu  dy 
( fo  P(u)  dv) 


(3.1)  dtu  +  .  V - (f  p' (u)  dtu  dy\  dyu  -  (  f  p(u)dy\  dyyu  =  0 

J0  p(u)  dy  \J o  /  \Jo  ) 

(3.2)  u(0,f)  =u(l,t)  =0, 

(3.3)  u(y,0)  =  u0(y)  =  9o(ys(0)). 


3.1.  Analysis  of  the  auxiliary  problem.  Let  a(u)  =  (  /0'  p(u)  dyj  and  let 

Ji(u)  be  the  solution  operator  for  the  linear  equation  dtw  —  a(u)dyyw  =  f.  More 
precisely,  w  =  W(u)/  implies 


dtw  -  a(u)dyyw  =  /, 
w(  0,  t)  =  io(l,  t )  =  0, 
w(y,  0)  =  w0. 


Now,  let  L  be  defined  by 


L(u)w  =  w  H — j- 


fo  P(u )  dV 


[  P‘(u)i 

Jo 


wdy  dyu. 
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L  is  a  rank-one  perturbation  of  I.  It  is  easy  to  check  directly  that 

_1  (fo  P'(u)s  dyj  y  dvu  (/q  p’(u)g  dyj  y  dyu 

’  ^  5  Jd  p(u)  d'y  +  fd  y  dyu  p’(u)  dy  9  p(u{l,t)) 

Therefore  by  (3.2)  L{u)  is  nonsingular  provided  p(u(l,t))  =  p( 0)  ^  0.  So,  if  u  is  the 
solution  of  (3.1)  then 


(3.4) 


u  =  T(u)  =  H(u) 


We  can  now  state  an  existence  theorem.  The  problem  is  solved  over  the  time  interval 
(0, T),  T  >  0.  As  usual,  for  a  nonintegral  positive  number  a,  C“([0, 1])  denotes  the 
space  of  Holder  continuous  functions  of  exponent  a  on  [0, 1].  The  corresponding  norm 
is 


l|(“)  — 
H(o,i)  - 


E  max  |  d 3xu(x)  |  +  sup 

j_0  x,x' €.(0,1), x^x1 


dxu(x)  -  dxu(x') 


where  [a]  stands  for  the  integral  part  of  a.  We  set  Q  =  (0, 1)  x  (0,T)  and  extend 
the  definitions  and  notations  to  functions  in  Q  in  the  obvious  way  [8].  A  general 
discussion  of  compatibility  conditions  of  the  type  used  in  the  next  result  can  be  found 
in  [8],  p.319.  Here  for  instance,  the  compatibility  conditions  of  order  1  for  (3.1),  (3.2), 
(3.3)  are 


«o(0)  =«o(l)  =  ug(0)  =<(1)  =0. 

Theorem  3.1.  Assume  the  density  function  p  is  of  class  C1 ,  is  positive,  uni¬ 
formly  bounded  and  bounded  away  from  zero.  Let  a  >  0  be  a  nonintegral  number. 
Then  ifu0  e  C3+“([0, 1])  and  satisfies  the  compatibility  conditions  of  order  +  1, 
the  problem  (3.1),  (3.2),  (3.3)  admits  a  unique  solution  in  C2+“,1+“//2([0, 1]  x  [0,T]) 
provided  \  p'  [oo  =  ma x.xeR  I  p'(x)  I  and  II  uo  1 1  (o~i)*^  are  sma ^  enoudh- 

Proof.  Let  Xo  =  C2+a,1+a/2(Q)  andA'i  =  C3+a,^~  (Q).  Using  classical  regularity 
results,  see  e.g.  [8],  p.320,  Theorem  5.2,  one  obtains  T  :  Xo  — t  X±,  where  T  in  defined 
in  (3.4).  Hence  T  is  a  compact  mapping  on  Xo .  Further 

||  T(u)  |Uc  <  C  (||  uo  ||SS“}  +  I  P'  looll  n  ||g+Q’1+“/2))  . 

Invoking  the  Schauder  Theorem,  we  get  existence  of  a  solution  in  X0  for  \p'  \oo  an(l 
||  uo  Hgt?  smaH  enough.  Uniqueness  is  easily  obtained  through  classical  arguments 
thanks  to  the  regularity  of  the  above  solution  via  an  adaptation  of  the  maximum 
principle,  see  [8],  p.22,  Theorem  2.8.  □ 

3.2.  Discretization  of  the  auxiliary  problem.  Let  Ax  =  1/(N  +  1)  and  let 

Xi  =  iAx ,  i  =  0, 1, . . .  JV+  1.  A  simple  finite  difference  approximation  of  (3.1)  is  given 
by 


N 


dtUi  +  4  I  E^m-W- 
\j=1 


2Ax 


(Ui+ 1  —  Pi- 1)  +  D(U)i  —  0,  i  —  1, . . . , N, 
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where  S  =  w jP{Uj ),  to/s,  j  =  0, 1, . . . ,  N  +  1  are  the  weights  of  the  numerical 

quadrature  and  where  Uo  =  f/jv+i  =  0.  Further,  in  the  previous  relation 

D(U)  =  ^  wjP(U j)  j  AU  e  Rn, 

where  A  =  {—dyy)h  €  RNxN  is  a  discretization  of  the  second  order  space  derivative 
operator  with  homogeneous  Dirichlet  boundary  conditions,  obtained  for  instance  by 
using  the  classical  three-point  formula,  and  U  is  the  N  x  1- vector  \Ui . . .  C/jy].  Let 
C  =  ( dy)h  €  RNxN  be  the  matrix  corresponding  to  the  discretization  of  the  convective 
term  in  (3.1),  here  Cij  =  (<5i,i+i  —  <5j+i)j)/(2Ax),  and  let  to  be  the  N  x  1-vector 
to  =  [wi . . .  wjv].  One  can  define  the  discrete  analog  to  L(u)  as 

Lh(U)  =1+  i(. xCU)(top'(U))T  e  RNxN. 

Note  the  outer  product  in  the  definition  of  Lh{U).  Using  the  new  variable  Z  =  dtU, 
i.e.,  Zi  =  dtUi,  i  =  1, . N,  the  semidiscretized  in  space  problem  can  be  written  as 
the  following  DAE 


(3.5) 

(3.6) 


dtU  =  Z 

0  =  Lh(U)Z  +  D(U). 


By  construction,  Lh(U)  is  equal  to  I  plus  a  rank-one  perturbation.  Elementary  Linear 
Algebra  implies 


LhiU)-1  =/- 


(. xCU){loP'{U))t 

S  +  (xCU)Tu>p'(U) ' 


In  other  words,  Lh(U)  is  nonsingular  provided  S+(xCU)T top1  (U)  7^  0.  As  was  the  case 
above  for  the  continuous  problem,  this  nonsingularity  condition  can  be  considerably 
simplified  and  turns  out  to  be  very  mild.  Indeed  after  summation  by  parts  and  some 
rearrangements,  one  obtains,  for  too  =  wjv+i  =  and  to±  =  . . .  =  tojy  =  Ax 


N+l 


N 


S  +  ( xCU)Ttop'(U )  =  £  tOjp(Uj)  +  • 


Uj+ 1  -  U 


3- 1 


j= 0 

1 


3= 1 


2Ax 


=  ^p(UN)+p(UN+1))  +  0(Ax\p" 


\su\2x), 


where  5U  is  the  vector  of  all  first  divided  differences.  Therefore,  the  corresponding 
nonsingularity  condition 


(3.7) 


\{p{UN)  +  P(UN+ 1))  +  0( Ax  |  p"  1 00  |  SU  t)  ±  0, 


is  clearly  a  discrete  analogue  to  p(0)  ^  0  since  f/jv+i  =  0. 

In  conclusion,  our  system  is  an  index-1  DAE  under  condition  (3.7).  We  apply  a 
linearly  Implicit  Euler  discretization,  see  e.g.  [5],  p.427  or  [2],  Chap.  4,  and  get 


(3.8) 


I 

-Atr 


-At  I 

-A  tLh(Un) 


jjn+l  _  jjn 

zn+ 1  -  Zn 


=  At 


Lh(Un)Zn  +  D{Un) 
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where  At  is  the  time  step.  To  find  the  expression  of  J" ,  we  need  to  take  the  derivatives 
with  respect  to  V  of  both  D(U)  and  Lfl(U)Z ,  since  J”  =  L'h(Un)Zn  +  D'(Un).  This 
can  be  done  by  direct  calculations;  those  expressions  are  omitted  here. 

The  matrix  to  be  inverted  in  (3.8)  is  nonsingular  for  At  small  enough.  For  the 
present  type  of  index- 1  DAEs,  the  usual  convergence  results  turn  out  to  be  still  sat¬ 
isfied.  One  easily  obtains  the  following  result. 

Theorem  3.2.  Let  ( U°,Z° )  be  consistent  initial  values,  i.e.,  Uf  =  9o(s(0)yi), 
i  =  1, . . . ,  N  and  Z°  =  —Lh(U°)-1D(U°).  Then,  under  condition  (3.7),  one  has 

II  U(tn)  -  Un  II  +  II  Z(tn)  -  Zn  II  =  O(At) 


Proof.  See  [3],  Theorem  1,  p.504.  Note  that  the  contractivity  condition  (1.4)  in 
[3]  is  trivially  satisfied  here  because  the  algebraic  equation  is  linear  in  the  algebraic 
variable  Z.  □ 

4.  The  full  problem.  Let  us  now  turn  to  the  powder  problem.  As  in  §3,  the 
problem  is  first  rewritten  in  a  constant  domain.  To  keep  the  notational  changes  to  a 
minimum,  the  original  non-transformed  variables  from  §2  are  denoted  with  a  bar,  and 
the  new  transformed  ones  are  taken  without  one.  The  relations  between  them  are 

y  =  z/H(t)  p(y,  t)  =  p(z,  t)  a(y,t)  =  a(z,t). 

When  rewriting  the  pressure  equation  (2.12)  in  terms  of  the  new  variables,  both 
H(t )  and  H'(t )  appear  in  the  terms  involving  time  derivatives.  Notice  that  M  = 
/(f W  l(d(z,  t ))  dz  =  fg  7 (a(y,  t))H(t)  dy,  and  thus 


(4.1)  H(t)  = 


M 


and  H'{t)  =  —M 


Jo  7 '{a)dtady 


Jo  7 {a(y,t))dy 
Having  “solved”  the  equation  for  H,  we  substitute  and  get 
(4.2) 


( fo  7 (°(y,t))dy) 


P 


7(c 

-dy  ( p( 


A  (7  (a)) 

A(y(a)dta) 

dtT(o)  +  — tt-t-ss  V  dyj(a) 


1 


7(a) 


b)i: 


a  ^  ,  A(i(a)dta)  _  , 

9tl{a)  +  A(7(a))  y9y 7(a)  J  dy 


~  "T^"^  (/  7^<T') dy )  dy^Kdyp^  =  °’ 

where  A{u)  =  -A  u(y)  dy.  The  stress  equation  becomes 

M 


(4.3) 


dva  +  dvp  + 


(— na  +  7)  =  0. 


fo  7(o-)  dy 

4.1.  Analysis  of  the  full  problem.  Define  the  integral  operators 


B(u)  = 


f  <z) 


dz  and  B(u)  = 


I"  u(z) 

J  0 


dz 
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The  boundary  conditions  on  the  transformed  variables  are 

9yp(0,  t)  =  0,  p(l,t)=patm  and  a(l,f)=0. 

The  initial  condition  on  p  is 

p{y,o)  =Po(y%  ye  (0,1). 

Hereafter,  the  above  problem  is  rewritten  as  a  parabolic  problem  for  the  stress  a,  and 
thus  an  initial  condition  ao  will  have  to  be  provided.  It  is  found  from  (4.3).  In  other 
words,  the  initial  stress  ao  is  taken  as  solution  to 

dya0  +  dypo  +  /(a0)  =0,  for  0  <  y  <  1, 

0o(l)  =  0, 

where  /(a)  =  fi  M  {l(°)  ~  k<j). 

J0  7 (<r(y))  dy 

Lemma  4.1.  Let  p0  e  C1([0, 1])  with  \p'0(y)  \  <  for  any  y,  0  <  y  <  1.  Then  if 
7  satisfies  (2.3),  there  exists  a  unique  function  a0  verifying  the  above  requirements. 
Proof.  We  construct  a  sequence  {crfi }„>0,  a°  =  0,  by  successively  solving 

dvaft  +  dyp0  +  fn{aft)  =0  for  0  <  y  <  1, 
ao  (1)  =  0, 

where  /„(a)  =  -jj — ■ — - (7(a)  ~  Ka)-  The  function  /„  is  Lipschitz  in  a  with 

J0  7«  (y))dy 

Lipschitz  constant  L  =  —  n).  Therefore,  there  exists  a  unique  function  aft, 

and  the  above  sequence  is  well  defined.  By  using  the  bounds 

,  .  kM  ,  ,  .  M  ,  .  „ ,  .  . 

/  (a)  = - —a  <  fn(a)  <  - (7  -  k a)  =  f+{a), 

7  7  m 

in  place  of  fn(a)  in  the  above  problem,  one  can  explicitly  construct  uniformly  bounded 
in  n  sub  and  supersolutions  to  the  above  problem.  Hence,  the  sequence  {o-q  }  is 
bounded  uniformly  in  n,  7  admitting  upper  and  lower  bounds  by  assumption,  see 
(2.3).  By  construction,  {dyaft}  is  consequently  also  uniformly  bounded.  Therefore, 
the  sequence  {aft }  is  equicontinuous  and  thus  by  Ascoli-Arzela’s  theorem,  one  can 
select  from  it  a  subsequence  converging  uniformly  to  a  Lipschitz  function  ao  which 
solves  the  above  problem. 

Uniqueness  follows  from  classical  arguments  due  to  the  local  Lipschitz  property 
of  /  as  a  function  of  a  and  to  the  regularity  and  boundedness  of  the  solution  to  the 
above  problem.  □ 

By  integrating  the  stress  equation  (4.3)  between  y  and  1,  0  <  y  <  1,  we  obtain 

(4-4)  p  =  P{a)  =  patm  -CT+  1  B{-na  +  7 (a)). 

A(7(a)) 

Taking  the  derivative  with  respect  to  t  yields 

(4.5)  dtp  =  -dta  -  ^^^("^  +  7(<0)  +  A(y(a))B^~Kdt<T  +  ^ ' 

The  next  step  in  the  analysis  consists  of  eliminating  one  of  the  two  unknowns  p  or  a 
to  obtain  an  integral  equation  similar  to  (3.1).  Here,  the  pressure  p  will  be  eliminated 
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to  obtain  a  problem  in  term  of  the  stress  a  only.  This  requires  some  manipulations 
of  integral  equations  involving  Volterra  and  rank-one  operators. 

By  a  Volterra  operator,  we  mean  an  integral  operator  of  the  form 

Vu(y)  =  [  v 
Jo 


(y,  z)u(z)  dz  or  Vu(y)  =  /  v(y,  z)u(z)  dz 


f 


where  the  kernel  v  is  assumed  to  be  in  L2.  It  is  well  known  that  the  spectrum  of  a 
Volterra  operator  is  the  single  point  0.  Rank-one  operators  can  be  expressed  as 


Ru(y )  =  (n  <g>  rr)u(y) 


rr(z)u(z)  dz, 


where  (-,  •)  denotes  the  L2( 0, 1)  inner  product. 

Lemma  4.2.  Let  V  be  a  Volterra  operator,  R  =  ri  <g>  rr  be  a  rank-one  operator, 
and  L  =  I  +  V  +  R.  Then  L  is  nonsingular  if  and  only  if 


(4.6) 


l  +  (rr,(/  +  V)-1ri)^0, 


in  which  case 


L-1 


(  {{I  +  V)-ln)®rr  \ 

V  1  +  (rr,(I  +  V)-Wl)J 


(i  +  v)-1. 


Proof.  Under  the  above  assumption,  (/  +  V)  is  clearly  nonsingular;  in  fact,  (/  + 
V)-1  is  equal  to  its  Neumann  series  1)"^".  Therefore,  the  result  is  a  simple 

consequence  of  the  Sherman-Morrison- Woodbury  formula,  see  e.g.  [4],  p.51.  □ 

Using  the  above  result,  one  can  thus  write 

(I  +  V  +  R)-1  -I  =  V  +  R, 


where  V  and  R  are  themselves  respectively  Volterra  and  rank-one  operators.  We 
define  the  following  Volterra  and  rank-one  maps 

v (u)  =  -  A(7(a))  B(~ku  +  T\v)u),  and  R(u)  =  +  7). 


By  setting  L  =  I  +  V  +  R,  (4.5)  takes  the  form  dtp  =  L(—dt<r).  The  rank-one  map  R 
can  be  written  as  R  =  ri  <g>  rr  with 


ri 


B(—k<j  -I-  7) 
M  A(,y(a))2 


and  rr  =  7  '(a). 


Since  by  (2.3),  the  kernel  of  V  is  bounded,  one  sees  that  the  nonsingularity  condition 
of  L  coming  from  (4.6)  reads  here  as  follows.  If  b  denotes 


(/  +  V)-in  =  £(-V)"n, 
n= 0 


i.e.,  the  unique  solution  of 
b- 


mwl (7>) 


—  n)b  dy  = 


A2  {n  (a)) 


[  (l(o-) 

Jy 


-  Ka)  dy, 
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then  L  is  not  singular  provided 

(4.7)  1  +  f  7'(ff)6dy^0. 

Jo 

The  next  step  consists  of  rewriting  the  full  problem  (4.2),  (4.3)  in  term  of  a  alone. 
The  terms  containing  unintegrated  time  derivatives  in  the  pressure  equation  (4.2)  are 

(1  -  7/T)%j  -  -  St7  =  -a(a)dta  -  (1  -  ^)  (V  +  R)dta. 

7  I 

where  a(a)  =  (1  —  7/r)  +  and  where  P  is  defined  in  (4.4).  Using  (4.3),  the 

diffusive  term  takes  the  form 

(1  -  ^)  A(7(a))2  dy(pKdyp)  =  -(1-1^1)  A(7(a))2  dy(P(a)Kdya)  +  a(a)  Q(a) 

with 

QH  =  (1  -  ^)dy(P(a)K( 7(a)  -  «a)). 

The  problem  can  then  be  rewritten  as 

a(a)L(dta)  =  (1  -  ^j^)  -4(7(a))2  dy{P{a)Kdya)  -  a(a)  Q(a), 


where  L  =  I  +  V  +  R  with 
1 


V(u)  = 


+ 


R(u)  = 


a(a) 

1 

a(a) 

1 

a(a) 


(1-3^)V(») 


l(c 


-  u)  I  B[i{cj)u  + 


ydyl(p) 


(1  -^W) 


+ 


1  P(q-)  A(y(a)u)  _ ]_  _  7(0).  ,  «  pr  x 

a(a)  7(a)  A(7(a))  ^  ofal  ^  r  ^  4^faU  ^  ^ 


a(a) 


^(7(0-)) 


We  have  implicitly  assumed  the  coefficient  a(a)  to  be  nonzero;  this  will  be  justified 
below  where  it  will  be  shown  under  what  condition  a(a)  is  positive.  The  rank-one 
map  R  can  be  expressed  as  R  =  ri  <3  rr,  where 


-  =  _J_ , ,  _  tW.  _J_  jV)  ydv7 _ (1  _  y  dvp 

1  a(a)  r  1  a{cr)  7(a)  ^4(7(0-))  a(cr)  T  A(^y(a))  ’ 

fr  =rr=  7' (a). 


Assuming 

(4.8)  l  +  t^a  +  y)-1^)  ^0, 

the  operator  Z  can  be  inverted  using  Lemma  4.2,  leading  to 


(4.9)  dta  -  rj(a)dy(P(a)Kdya)  +  Q(a)  =  F(a), 


POWDER  CONSOLIDATION 


13 


where 

V(v)  =  -7-t  (1  -  -^r-)  ^{o))2 , 

a{a)  1 

F(a)  =  (V  +  R)  (r,(a)dy(P(a)Kdya)  -  Q(a)  j  . 

In  the  above  equation  (4.9)  for  a,  it  is  easily  checked  that  the  diffusive  coefficient  is 
positive  provided  P(a)  stays  positive.  However,  as  can  be  seen  from  (4.4),  P(c r)  is 
positive  only  for  small  enough  values  of  a,  due  to  the  growth  condition  on  7,  see  (2.2), 
(2.3).  Accordingly,  we  first  regularize  (4.9)  to  ensure  parabolicity,  establish  existence 
of  a  solution  to  the  corresponding  problem,  and  show  under  what  conditions  the  result 
extends  to  the  original  nonregularized  problem. 

Let  (0 ,  T),  T  >  0,  the  interval  in  which  the  problem  is  to  be  solved.  Considering 
P  from  (4.4)  as  an  operator  from  C([0, 1]  x  [0,  T])  into  itself,  we  define  a  “regularized” 
Pe  as  follows.  If  a  stands  for  the  unique  solution  of  7 (a)  =  na,  see  Figure  2.1,  then 
we  set  <7*  =  min (p0tm,d).  Further,  for  any  a  £  C{[ 0, 1]  x  [0,T])  and  for  some  e  >  0, 
we  define 


St  =  {{y,t)  €  [0, 1]  x  [0 ,T];a(y,t)  <  a*  -  e}. 

The  regularized  function  Pe  is  taken  as 

Pt(a)\S:  •  =  P(a)  ,s> 

Pe(a)>e,  Va  S  C([0, 1]  x  [0,  T]), 

Pe  is  a  smooth  function  of  a. 

The  corresponding  regularized  problem  is  then 

(4.10)  dta  -  r)PeKdyya  -  Vdy(PeK)dya  +  Qe(a) 

(4.11)  dya{  0,  t)  +  — (j(cr(0,  t ))  -  k<t(  0,  t)\ 

Jo  7Wy,t))dy\  J 

(4.12)  a(l,t) 

(4.13)  a(y,  0) 

with  Qe  and  Fe  are  defined  are  in  a  natural  way,  P  being  replaced  by  Pe.  The  condition 
(4.11)  at  y  =  0  is  a  direct  consequence  of  (4.3)  and  dyp(0,t)  =  0. 

Lemma  4.3.  Let  0  <  a  <  1,  and  let  the  nonsingularity  conditions  (4-7)  and  (4-8) 
be  satisfied.  Let  also  the  assumptions  of  Lemma  4-1  be  verified,  assuming  further  that 
p0  £  C3+“([0, 1]).  Then,  if  the  initial  stress  op  €  C3+“([0, 1])  from  Lemma  4-1  satisfies 
the  compatibility  condition  of  order  1,  the  regularized  problem  (4-10)-(4-13)  admits  a 
solution  <7e  in  C2+a,1+“/2([o,  1]  x  [0,T])  provided  7^,  Ko,  T  and  ||  op  are  small 

enough. 

Proof.  As  was  done  for  the  auxiliary  problem,  a  linearized  solution  operator 
corresponding  to  the  above  regularized  problem  is  constructed.  More  precisely,  for  a 
given  u  £  X0  =  C2+“’1+“/2(<2),  consider  Te(u)  =  a  where  a  satisfies 

dtcr  -  r](u) P€(u) K (u)dyya  -  r](u)dy(Pe(u)K(u))dyu  +  Qe{u)  =  Fe{u), 

<V( 0,  t)  +  M  (iMO,  t))  ~  kct{ 0,  t)\  =  0, 

Jo  T{u{y,t))dy  \  J 


=  W, 
=  0, 

=  0, 

=  °o(y), 
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A  dyo(\.  t )  +  ct(1,  t)  —  A  dyu{\,  t)  =  0, 

a(y,0)  =  a0{y), 


In  order  to  avoid  having  to  treat  a  mixed  Dirichlet-Robin  problem,  the  given  condition 
at  y  =  1,  i.e.,  a(l ,t)  =  0,  has  been  rewritten;  the  coefficient  A  ^  0  is  to  be  chosen 
below.  Note  that  when  u  =  a  is  a  fixed  point  of  Te,  then  a(l ,  f)  =  0  and  thus  the 
boundary  value  is  correctly  recovered. 

Using  classical  regularity  results,  see  e.g.  [8],  p.  320,  Theorem  5.3,  one  obtains 
the  existence  of  a  unique  solution  a  =  Te(u)  to  the  above  linearized  problem.  Further 

II  a  |Uc  =  ||  Z(u)  |Uc  <  c(||  a0  ||gj“}  +  K0  (||  u  ||g+«,1+«/2)) 2 

+  ^(F  +  i'oo)  (ll  «(0,  •)  ||g+f 72  +  ||  u(  1,  •)  ||g+“)/2)) , 

where  A  was  chosen  so  that  the  contributions  from  both  boundary  terms  “match”  in 
the  above  inequality.  Therefore,  Te  maps  Xo  in  X\  =  C3+“,_s_  (Q),  and  is  thus  a 
compact  mapping  on  Xo .  The  existence  of  a  fixed  point  ae  =Te(&e)  results  from  the 
Schauder  Theorem.  The  function  ae  clearly  satisfies  (4.10)-(4.13).  □ 

Note  that  the  restriction  a  <  1  in  the  statement  of  Lemma  4.3  is  not  essential. 
Indeed,  higher  values  of  a  can  be  considered  as  the  price  of  having  to  bound  higher 
derivatives  of  the  function  7  and  satisfying  a  compatibility  condition  of  order  1  +  [f  ]  • 

Clearly  now,  using  again  the  classical  regularity  result  invoked  in  the  latter  proof, 
for  small  enough  data,  the  largest  value  of  ae  will  not  exceed  a*  —  e,  threshold  above 
which  the  regularization  of  P  takes  place.  Uniqueness  of  the  solution  follows  from 
the  same  argument  as  in  the  proof  of  Theorem  3.1.  We  have  thus  established  the 
following  existence  result  for  the  main  problem. 

Theorem  4.4.  IfT  >  is  small  enough  and  under  the  assumptions  of  Lemma  f.3, 
there  exists  a  €  C2+“,1+“//2(H),  p  €  C2+“’1+a//2(H)  and  H  €  C1+a//2([0,  T])  which 
satisfy  (2.12)-(2.16),  where  S  =  {(z, f);0  <  t  <  T,0  <  z  <  H(t)}. 

4.2.  Discretization  of  the  full  problem.  The  problem  is  discretized  under 
its  formulation  (4.2),  (4.3),  i.e.,  we  work  with  the  transformed  variables  in  a  fixed 
rectangular  domain  (0, 1)  x  (0,T).  Due  to  the  different  type  of  boundary  conditions, 
the  mesh  related  quantities  are  slightly  modified  from  §3.  We  set  Ax  =  1/N  and 
let  Xi  =  (i  —  1) Ax,  i  =  1 ,N  +  1.  The  semidiscretized  variables  are  the  pressure 
P(t)  =  [Pi(t), . . . ,  Pjv(t)]  and  the  stress  S(f)  =  [Ei(f), . . . ,  Sjv(i)]-  The  boundary 
conditions  at  y  =  1  read  P/v+i  =  Patm  and  Ejv+i  =  0.  We  introduce  the  following 
discretized  operators 

Cp :  N  xN  matrix  corresponding  to  a  second  order  centered  discretization  of  dy  with 
pressure  boundary  conditions  (at  y  =  0  and  y  =  1), 

D:  N  x  N  matrix  corresponding  to  a  second  order  centered  discretization  of  dyy 
with  pressure  boundary  conditions  (at  y  =  0  and  y  =  1), 

Ca :  NxN  matrix  corresponding  to  a  second  order  essentially  centered  discretization 
of  dy  with  stress  boundary  condition  (at  y  =  1). 
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The  construction  of  Cp  and  D  is  elementary.  For  Ca.  we  take 


C„ 


1 

2Ax 


-3  4  -1 

-1  0  1 
0-1  0 


0 

0 

1 


0 

0 

0 


0  ...  0  -1  0  1 

0  ...  01/6  -4/3  1/6 


We  also  introduce  discretized  counterparts  to  the  integral  operators  A  and  D 


IV+l  i 

-4a W  =  -  E  u>wi  and  (^A (W))i  =  E  UjWj  • 

1  j= i  j= i 

Finally,  we  introduce  some  notation  for  the  time  derivative  of  the  two  main  unknowns 
P  and  E  by  setting  U  =  dtP  and  V  =  dp.  The  discretized  problem  is  then 


(4.14)  dtP  =  U 

(4.15)  o  =  (i-I)(c/  +  MA)rcpp|- 


-P~-\)crp-?Ac^ 


6MV)  +  ^AjA1ba(~i'yc,s) 


+(1  -  ^Ma(7)2£CP,  S)P  +  bcl  =  f(P,  U,  V,  E) 


(4.16)  0  —  -|-  CpP  H~  bc2  + 

(4.17)  ftE  =  y, 


-4a  (l) 


-4  a  (7) 
S) 

(-«E+7)=5(P,C/,y,S) 


where  7  and  7'  are  to  be  understood  as  the  vectors  7(E)  and  7'(S)  and  where 
Y  =  [xi, . . .  The  two  vectors  be  1  and  bc2  result  from  the  presence  of  bound¬ 
ary  conditions.  Finally,  the  vector- vector  multiplications  in  the  above  expressions  are 
to  be  understood  component  by  component. 

The  structure  of  (4.14)-(4.17)  is  slightly  more  complicated  than  that  of  the 
semidiscretized  problem  (3.5),  (3.6).  More  precisely,  (4.14)-(4.17)  has  the  form 


dtP  =  U, 

0  =  f(P,  U,  V,  E):., 
0  =  g(P,  0, 0,  S), 
dp  =  V. 


The  above  system  corresponds  to  a  semi-explicit  index  2  DAE  or  equivalently  to  a 
fully  implicit  index  1  DAE  [5],  § VII.3,  VII.4.  Indeed,  it  is  clearly  equivalent  to 

(4.18)  0  =  f(P,  dtP,  dp,  S), 

(4.19)  0  =  g(P,  0, 0,  S), 

which  is  the  discrete  counterpart  to  (4.2),  (4.3).  Note  that  gi  =  Cp  being  nonsin¬ 
gular,  one  can  solve  (4.19)  into  P  =  P(S),  which  corresponds  to  (4.4).  Now,  by 
differentiating  (4.19),  we  get 


gidtP  +  gtdp  —  0. 
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This  leads  to 


(4.20)  dtP  =  -CplgAdtX, 

which  is  a  discretized  version  of  (4.5).  Using  again  the  Sherman-Morrison- Woodbury 
formula,  one  can  check  that  gA  is  also  nonsingular  provided  that 


(4.21) 

ca  + 

(4.22) 

1  -  - 
Tl 

1 


-4  a  (7) 

1 

MAA(-f)2 


—kI  +  diag(^'  (Y))j  is  not  singular, 

(7(E)  -  «S)T(^7'(S))  A  0. 


It  is  easy  to  verify  directly  from  the  matrix  structures  that  for  Ax  small  enough  (4.21) 
is  satisfied.  Condition  (4.22)  is  the  counterpart  here  to  (3.7)  for  the  auxiliary  problem 
and  (4.7)  for  the  nondiscretized  full  problem.  Note  that  it  is  satisfied  if  max*  ^77^  is 
small  enough.  Finally,  plugging  (4.20)  into  (4.18),  we  get 


0  =  /(P(S),-C'-V^S,^S,S). 


Under  an  involved  nonsingularity  condition  that  is  the  discrete  counterpart  to  (4.8), 
the  latter  equation  can  be  solved  for  dt S.  It  is  thus  a  fully  implicit  index  1  DAE. 
Results  similar  to  the  one  use  in  §3.2  to  establish  error  estimates  can  also  be  found 
for  the  present  type  of  problems,  see  e.g.  [11],  [5]  Chap.  VII.  We  do  not  pursue  this 
issue  further  in  this  paper. 

The  time  discretization  of  (4.14)-(4.17)  is  again  taken  as  linearly  Implicit  Euler, 
which  reads  here 


I 

-At  I 

0 

0 

pn-\- 1  _  pn 

-  jjn  - 

-At  A” 

-At /2” 

-At /3" 

-At  A” 

un+ 1  -  un 

=  At 

r 

-At  g™ 

-At  A1 

-At  3? 

-At  A1 

yn+i  _  yn 

9n 

0 

0 

-At  I 

/ 

^n+1 _ 

yn 

In  the  previous  relation,  the  subscripts  denote  derivatives  with  respect  to  the  corre¬ 
sponding  variables.  All  the  first  partial  derivatives  of  /  and  g  are  needed.  For  /,  we 
note  that  /"  =  diag(  1  —  7(S")/r);  the  rest  of  the  terms  are  evaluated  numerically 
through  finite  differences.  For  g,  we  obtain 

9i=CP 

92  =  9s  =  0 

3?  =  Ca  +  (-«/  +  diag( 7'(S")))  -  \  (-kZ  +  7 (S))(W7'(S))T. 

-4  a  (7)  MAa(7)2 

Note  the  outer  product  in  g " .  The  position  of  the  free  boundary  can  be  determined 
a  posteriori  through  (4.1),  i.e.,  Hn  = 

5.  Computational  results.  We  now  illustrate  the  feasibility  and  efficiency  of 
the  above  numerical  method  by  discussing  some  computational  results.  The  bulk 
density-stress  relationship  is  taken  as  (2.2).  The  following  values  of  the  various  pa- 
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pressure  (lbf/ft2) 

1 2210 
1 2200 
L  1 2190 
-  -  2180 

-  -  2170 

-  -  2160 

-  -  2150 


0  100  200  300  400  500 

time  (s) 


stress  (lbf/ft2) 


time  (s) 


Fig.  5.1.  Calculated  pressure  and  stress  fields  in  the  original  geometry;  N  =  100,  NT  =  100, 
T  =  500s. 


rameters  have  been  used  and  were  chosen  so  as  to  correspond  to  a  realistic  situation 

P  m  =  -25 
7 m  =  60  lbs/ft3 

am  =  13  lbs/ft2 
7o  =  80  lbs/ft3 
T  =  200  lbs /ft3 
K0  =  10~4  ft4lbs  xs  1 
a  =  4 
k  —  1  ft-1 ; 

the  atmospheric  pressure  is  patm  =  2116.2  lbs/ft2.  The  problem  is  initialized  as 
follows.  The  initial  height  of  the  powder  column  is  prescribed;  here  H( 0)  =5  ft.  The 
bunker  is  assumed  to  have  a  unit  cross  section  (i?  «  .56ft).  A  initial  pressure  field  p0 
which  is  consistent  with  the  boundary  conditions  in  (2.16),  is  given  as 

B,W=  (1  +  lTO  1 3  ~  (ff(o) ,2))  ' 

In  the  calculations  below,  we  use  N  =  100  and  NT  =  T/At=  100. 

Figure  5.1  gives  a  typical  illustration  of  the  behavior  of  the  problem.  Note  that 
Figure  5.1  is  related  to  the  original  geometry  of  the  problem,  i.e.,  the  variable  is  z 
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stress  (lbf/ft2) 


0  500  1000  1500  2000 


time  (s) 


height  (ft) 


Fig.  5.2.  Calculated  stress  field  in  the  original  geometry;  left:  calculation  from  0  to  T  =  2000s, 
right:  asymptotic  stress  profile;  N  =  100,  NT  =  400. 


not  y.  As  expected,  one  observes  a  settlement  of  the  powder.  For  both  graphs,  the 
top  line  corresponds  to  the  position  of  the  top  of  the  powder  column,  i.e.,  the  free 
boundary.  Further,  the  pressure  field  is  found  to  asymptotically  converge  to  a  uniform 
pressure  value  corresponding  to  patm ,  the  atmospheric  pressure,  i.e.,  equilibrium  of 
the  pressure  is  established.  The  stress  field  is  also  found  to  converge  to  a  stationary 
distribution  a which,  in  term  of  the  original  variable  z,  is  solution  to 

A  (  rcooo  Tq^croo))  —  0,  0  z  <c  c, 

^00(^00)  =  0; 

where  H ^  is  the  asymptotic  value  of  the  height  of  the  powder  column.  From  Fig¬ 
ure  5.1,  one  can  observe  that  the  powder  has  not  fully  settled  after  500s.  Figure  5.2 
illustrates  the  convergence  to  the  asymptotic  stationary  state  after  2000s. 

6.  Conclusion.  A  mathematical  model  of  the  phenomenon  of  powder  consolida¬ 
tion  has  been  derived  and  analyzed.  A  robust  numerical  method  has  been  proposed, 
and  successfully  implemented.  Several  questions  and  generalizations  deserve  further 
study. 

The  model  is  a  generalization  of  Janssen’s  approach,  but  it  does  rely  on  the  same 
two  crucial  assumptions  (quasiuniformity  on  horizontal  cross  sections,  and  horizon- 
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tal/ vertical  alignment  of  the  principal  stresses)  that  may  be  questionable  in  some 
cases.  How  to  bypass  those  assumptions  will  be  the  object  of  future  work.  From  a 
more  applied  view  point,  the  bunker  containing  the  powder  may  be  axisymmetric  but 
of  variable  cross  section.  This  aspect  may  be  easily  included  in  the  present  model. 
Further,  one  may  consider  adding  some  powder  from  the  top  and/or  retrieving  some, 
typically  through  outlets  at  the  bottom  of  the  bunker.  While  the  first  case  seems  to 
bring  in  some  modeling  difficulties,  the  second  one  (retrieval)  just  amounts  to  chang¬ 
ing  one  boundary  condition  and  can  thus  be  handled  without  difficulty  in  the  present 
approach. 
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